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CHAPTER i 


INTRODUCTION 

1.1 INTRODUCTION TO DECONVOLUTION : 

The output of a linear system can be represented by the 
convolution of the input and the system impulse response. The system 
impulse response itself can, in many cases, be represented as a 

convolution of impulse response functions of many subsystems. 

Deconvolution is the process of removal of any or all of the 
convolutional factors. 

Consider a linear system with impulse response , h(n). The output 
of the system to an input x(n) is given by 

y(n)» x(n)* h(n) 

There are many applications in which given y(n) and either x(n) or h(n) , 
we wish to determine the unknown signal. In deconvolution process , we 
wish to recover either h(n) or x(n) from the given information. 

Deconvolution plays an important role in identification of physical 
systems. We excite the system with a known input , then measure the 

output and obtain h(n) from a knowledge of x(n) and y(n). Let q(n) be a 

function such that 

x(n)*q(n)= S(n) it.l) 

where $( n) is the unit impulse sequence. Then convolving y(n) with a(n) : 
we get 

y(n)f q(n) = h(n) * x(n) * q(n) = h(n) * S(n) * h(n) 

The function q(n) is called the deconvolution filter , sinoe it ; 



"deconvolves” the input x(n) from the output y(n) 


Another variation of the deconvolution problem is to find a linear 
system with impulse response g(n) such that 

y(n)* g(n) = x(n) 

which implies that 

h(n)*g(n) = S(n) (i.2) 

The system g(n) is known as the "inverse" of the system h(n). The 
inverse systems find applications in situations where it is necessary to 

compensate for the (undesirable) effects of a linear system. This is 

commonly referred to as "channel equalization". 

The problem of deconvolution is concerned with the development of 
fast and stable algorithms to implement deconvolution / inverse filters . 

Deconvolution finds applications in diverse areas. Deconvolution of 
speech can provide the useful information of the vocal tract. In seismic 
signal processing , the deconvolution of a seismogram can be used to 
estimate the characteristics of the earth's subsurface. In passive sonar 
signal processing, deconvolution of the received signal can be used to 
estimate the target's radiated signal. 

12 METHODS OF DECONVOLUTION : 

There are a number of deconvolution schemes that are currently 
being used. Each method is based on certain models and assumptions. So 
the success of a deconvolution method depends upon how well the 

assumptions made by it are met in the application. There are three 

popular methods, namely, Predictive Deconvolution, Homomorphic Deconvolution 
and Deconvolution by Kalman filtering, which are briefly outlined below . 



i. 2.1 


PREDICTIVE DECONVOLUTION Cl], LSI, CIO], C15], C17], [21] 


This makes use of the following statistical model . The system is 
characterized by a random impulse response h(n) given by , 

h(n) = S a i n - m i ) (1.3) 

i 

where is a randomly occuring point in time that can take on discrete 
values CO, 1, } ; and the a^'s are uncorrelated. 

A minimum phase signal x(n) is assumed to be the input to the 
above system. Under the uncorrelated random feature of the system 
impulse response, the autocorrelation, 4>yy , of the output y<n) is given 
by , 

fyy ( n) = * xx (n) * f> hh <n) 

= #xx<n) * K.S(O) , K is a constant. 

= K.* xx (n> (1.4) 

Thus the autocorrelation, 4>yy , of the output y(n) is equal to the 
autocorrelation # x x<n) of x(n) multiplied by a constant. This scaling has to 
effect on deconvolution procedure and hence can be ignored. The 4>yy 
can be computed directly from the output y(n), and thus <P XX can be 
obtained. The autocorrelation of x(n) is sufficient to determine the 
deconvolution filter for x(n) , viz., the prediction error operator with 

prediction distance equal to unity. This prediction-error operator, when 

convolved with x(n>, produces a spike at n=0, and hence is called the 

deconvolution-operator. Then the system impulse response h(n) can be 
computed by convolving this prediction error operator with the output 
y(n). The prediction filter with prediction distance equal to unity is 



called Wiener -filter . 


i. 2. 2 HOMOMORPHIC DECOW VOLUTION til, [33, [53, [133, [233 : 

Homomorphic systems are a class of nonlinear systems which satisfy 
a generalized principle of superposition. Such systems are particularly 
useful in separating the signals combined through convolution. Consider 
the transformation defined by, 

y = T (x). 

If T is a linear system, it satisfies the superposition principle 

defined by , 

T ( ax l + b %2 ) = a T(Xj() + b T(x 2 ) (1.5) 

In Homomorphic signal processing, a system D, referred to as the 
characteristic system, transforms a convolutional space to an additive 

space. It is defined by the relation, 

D ( axj * bx 2 ) = a D(xj) + b D(x 2 ). (1.6) 

The system, D"^, the inverse of D performs the transformation from 


an additive 

space back 

to 

the 

convolutional space. 

If in the additive 

space, the 

contribution 

due 

to 

the two convolutional 

factors are well 


separated , one of the components can be removed by low-pass or high- 
pass filtering in the additive space ( or cepstral domain ). 

1.2.3 DECONVOLUTION BY KALMAN FILTERING £13 

The deconvolution problem can be solved either by inverse filtering 
or calculating the prediction error. Both the methods are well known in 
the theory of Wiener filters. If, however, the generating process of the 
signal is known and can be described by a set of differential equations, 



then Kalman filter can also be used to solve the deconvolution problem. 
The modeling of the linear system in Kalman's formulation is based on a 
state space representation of linear time varying systems via a first 
order vector differential equation instead of a time-varying impulse 
response, as in the case of Wiener filter. This formalism results 
in specification of optimum estimates as the solution to a differential 
equation whose coefficients are determined by the statistics of the 

processes involved. 

Kalman filtering approach to deconvolution is applicable to time- 
varying or time-invarient signal as well as to nonstationary or stationary 
noise processes. One interesting feature of this approach is the 

reversal of the intuitive role of the input and the system for the same 
model as assumed by Predictive Deconvolution. In this case the signal x(n) 
is considered as the impulse response of the system and the random 
sequence, h(n), is considered as the input and is modeled as noise. 

1. 2. 4 OTHER METHODS : 

Besides the three techniques discussed above, there are many other 

methods used for deconvolution. The scheme called, norm method of 

deconvolution is suggested by Claerbout and Muir £71. This method is said 
to have robustness characteristics. The Minimum Entropy method of 
deconvolution suggested by Wiggins C253 is similar to Predictive 
Deconvolution in the sense that operators are determined direotly from 
the data, but Predictive Deconvolution seeks to maximize entropy as 
opposed to the former method, which minimizes the entropy. In the 
Deterministic Deconvolution method £1] , the input to the system is given 



and is deconvolved from the output of the system. A time-varying 
deconvolution method has been developed which is based upon adaptive 
linear filtering techniques. In this Adaptive Deconvolution scheme Cil , 
filter coefficients are designed for each sample of the input trace using 
an adaptive algorithm. The Dynamic Predictive Deconvolution method C183 is 
used for the system which is modeled as a lattice filter. 

i.3 RELATIVE MERITS OF THE DECONVOLUTION METHODS: 

The Predictive Deconvolution is very effective when the signal to 
be deconvolved is minimum phase . If the signal is nonminimum phase, then 
Homomorphic filtering or Kalman filtering can be applied. Further, 
Homomorphic Deconvolution does not require one of the convolutional 
factors to be random in nature, whereas for Predictive Deconvolution it 
is an essential requirement. However, the success of Homomorphic 
Deconvolution depends upon the degree of separability of the convolutional 
factors in the cepstral domain. Kalman filtering can be used in time- 
varying cases. 

i . 4 APPLICATIONS OF DECONVOLUTION 
1. 4. 1 SEISMIC SIGNAL PROCESSING : 

Deconvolution is applied in seismic signal processing to estimate 
the physical properties of the earth's subsurface, from which the 
presence of oil and natural gas in the earth's subsurface can be 
predicted. A major goal in the seismic deconvolution is to remove all 
the reverberation effects to retrieve the earth's reflectivity function. : 


I 


1 



1. 4. 2 SPEECH PROCESSING 


Deconvolution is used in speech analysis, e.g. in an automatic speech 
recognition system, we begin with the speech waveform and the desired 
result is the recognition of the speech. Other applications are speaker 
identification and speaker verification. Another application is secure voice 
transmission and data rate compression, which involves speech analysis 
followed by speech synthesis. If the speech is transmitted in PCM form , 
the data rate required of the order of 90000 bits per second. Through 
the use of speech analysis using deconvolution, followed by appropriate 
coding, transmission and resynthesis at the receiver this can be reduced 
by a factor between 10 and 50. 

1. 4. 3 SONAR SIGNAL PROCESSING : 

In passive sonar operations, a target radiates an acoustic signal 

into an accoustic medium with surface and bottom boundaries. Generally, 
the boundaries of this medium produce multipath interference. 

Consequently, a passive receiver receives a distorted version of the 
radiated signal. Deconvolution can be used to suppress the effect of 

multipath interference in the received signal. 

i. 4. 4 OTHER APPLICATIONS 

Deconvolution has a number of other applications. Examples include 

the probing of dielectric materials by electromagnetic waves, the study of 
optical parameters of thin films, the probing of tissues by ultrasound 

and the design of broadband termination of transmission lines. 


1.5 


MODELING IN SPEECH AND 


SEISMIC SIGNAL PROCESSING 

Every deconvolution scheme is based on a model of the generating 
process. A model is a representation of a process in which the details 
that appear unessential for the intended use are omitted. In speech 
analysis, the speech waveform is modeled as the response of a linear 
time- varying system, the vocal tract, with the appropriate excitation. In 
seismic signal processing, the seismogram is modeled as the response of 
a linear time-invarient system, the earth, to a seismic wavelet as the 

input. So in this chapter, the modeling of the vocal tract and the 
earth as a " lattice filter " is considered. 

1. 5. 1 MODELING OF SPEECH WAVEFORM : 

Speech is produced by excitation of the vocal tract by glottal 

pulses C1S1. If the vocal tract shape is fixed, it is modeled as a linear 
time-invariant system and the output is the convolution of the excitation 
and the vocal tract impulse response. Different sounds are produced by 
changing the shape of the vocal tract. If the vocal tract shape changes 
slowly, for short periods of time ( of the order of 20 - 30 milliseconds), 
it may be assumed to maintain a fixed configuration. 

In speech, the vocal tract is modeled as an acoustic tube of 

varying cross sectional area. It is considered to to be a set of 

interconnected acoustic tube sections of equal length and of different 
cross sectional area < Fig. i.i ). It is generally assumed that sound 
propagation through each section can be treated as a plane wave, and 
the internal losses can be ignored. The acoustic impedance of sound 


wave varies inversely with the tube area, i.e. Z = ( p c) / A , where, fi, c, 
A are the air density , the velocity of the sound and the tube area 
respectively. 

Therefore, as the sound wave propagates from the glottis to the 
lips, it will suffer reflections every time it encounters an interface ; 
that is, every time it enters a tube segment of different diameter ( Fig. 
1.2). So within each section, there will be forward and reverse travelling 
waves ; the forward travelling wave moves in the direction from glottis 

to the lips , and the reverse travelling wave moves in the direction 
from lips to the glottis. At the boundary between each section , some 
fraction of the forward travelling wave gets transmitted to the next 
section , and some fraction is reflected back, the same is true for the 
reverse travelling wave in each section. Multiple reflections will be set 

up in each segment and the tube will reverberate in a complicated 
manner depending on the number of segments and the diameter of each 

segment . 

Let i± m be the reflection coefficient at the boundary between mth 

and ( m—l)th sections and u + m and u~ m be the forward and reverse 

travelling waves respectively , & " t " denote half the time required for 

a wave to travel from one end of the section to the other. It can be 

shown that 

u+ ro-l < 1 + r > = Vm u ~ m -t < t - t ) + < i + jti m ) u + m C t - t ) (i.7) 

and 

u“ m ( t + 7 ) = ( i - ) u-^j < t - t ) - ti m u+ m ( t - t ) 

The above two relationships are described by a lattice filter , with the 
lattice parameters being equal to /ij ' s , the reflection coefficients of 



the acoustic tube model of the vocal tract ( Fig. 13 ). For this lattice 
filter structure , the delay time 2r of each section is taken to be a 
unit delay. As the shape of the vocal tract changes during utterances 

of different sounds , the cross sectional area of each section of the 
model , or equivalently the reflection coefficients (Mm , are modified. 

The speech waveform coming out of the tube is recorded and the 
tube parameters such as the reflection coefficients and the width of 

each section are extracted after deconvolution. 

1. 5. 2 MODELING OF SEISMOGRAM [23, C93, C173, C201, [213: 

The mathematical model for the real world seismic deconvolution 

problem is three - dimensional. This three - dimensional subsurface model 
is quite difficult to analyze ; it is still in research stage. So a one - 
dimensional subsurface model is considered. First , the fundamentals of 

seismic data processing are described. Then the modeling of the 
earth is discussed. 

1.5.2. 1 Seismic data processing : 

Deposits of petroleum and natural gas are located in reservoirs 

deep below the surface of the earth. In order to find these resources 

geophysical exploration methods are used. Exploration seismology may be 

broadly divided into the reflection seismology and the refraction 

seismology. Most of the exploration of petroleum and natural gas is 
done by reflection seismic methods. 

Reflection Seismology : 


A source at the surface is set off. The source may be a 


dynamite explosion or dropping of a weight , or some other means to 
transmit seismic energy into the earth. The seismic waves travel from 

the source to the subsurface layers within the earth where they are 
reflected. The reflected seismic waves are then recorded by instruments 

on or just below the surface. This recorded signal is known as a 
seismic trace or seismogram. The earth's sedimentary layers are 

approximately horizontal , but they do have features such as anticlines 
that can serve as traps for petroleum and natural gas. 

Each seismic trace is a time series made up of reflected waves 
together with various interfering waves and noise (Fig 1.4). The desired 
events are the primary reflections and the undesirable interfering waves 
are due to the multiple reflections. In order to suppress this 

interferance, redundancy in data acquisition is incorporatd. So in reflection 
seismology , seismic data are collected in a special way. A single impulse 

source is used along with many receivers, which are spaced eaually 
along a line, for recording the seismogram. The source is activated and 

the traces are recorded. Then the entire configuration of the source 


and the 

receivers is 

moved by 

a distance equal 

to 

one half the 

distance 

between two 

receivers , 

and another set 

of 

seismograms are 


recorded for second source excitation. This procedure is repeated a 
number of times. By moving the configuration in such increments, each 
depth point is covered several times , and thus a multiple coverage of 
a depth point is obtained. Reoording of seismic data by a multiple 
coverage scheme introduces considerable amount of overlap or redundancy. 
«AU the traces thus recorded , can be sorted out ( or gathered ) into a 
group such that all traces within that group have a common depth point 



This is known as CDP gathering (Fig i.5 ). 


Next , some static and dynamic corrections are made to the seismic 
traces , in each group. A common depth - point gather or group is taken. 
Because of lateral variations in the thickness of the near - surface 
layers , each trace is corrected by a time shift refered to as static 
correction. The static correction has the effect of placing source and 
the receiver on a fictitious horizontal datum plane. The dynamic correction 
is to overcome the effect of offset between source and receiver . 
Increasing the offset increases the arrival time of a reflection from a 
given interface as a result of increased distance travelled by the wave . 
So the dynamic correction is the correction for the time of first 
arrival of the wave arriving at each receiver so as to convert each 
trace to the equivalent trace that would have been received if both the 
source and the receiver are at the same point ( i.e. at a point directly 
above the common depth point of the group) . In other words, under the 
appropriate dynamic corrections , the primary reflections of all the traces 
in the group will be in phase, thus making the corrected traces 


coherent. Then next 

step is to 

add 

all the 

traces in 

a group. By 

adding, the primary 

reflections , 

which 

are in 

phase they 

are amplified. 


whereas the random noise and the multiple reflections are suppressed. 
The resulting trace is equivalent to the one , which could have been 
obtained if the source and the receiver were at the same point, and is 
called " normal incidence seismogram " . 

After the signal enhancement , the next step is to reoover the 
information about the earth's subsurface. This can be done by 

deconvolution of the seismogram ( i.e. the normal incidence seismogram ). 


i. 5.2.2 Modeling of the earth 


There are two basic modeling approaches in seismic data 
processing. 

1. Random fflfldel : 

In the Random model , the depths and the reflectivities of the 
deep reflecting horizons are considered to be random in nature. This 

consideration is valid in seismology because, the widths of the different 

layers and their composition are different. As the reflection coefficient at 

the interface depends upon the the acoustic impedance , which inturn 
depends upon the composition of two adjoining layers , the reflection 
coefficients can be considered to be random. The source signal is not a 

unit spike but a minimum or mixed delay signal called seismic wavelet. 

The surface layers produce reverberation due to multiple reflections. 
Thus the seismogram is the convolution of the seismic wavelet , the 
random reflection coefficient series and the reverberation sequence. 

2. Dynamic model : 

The earth's geological layers are characterised as elements of a 

stratified system for the purpose of modeling the acoustic wave 
propagation in various layers. The dynamic model of the earth is based 
on the following facts : wave motion in each layer is characterised by 
two components travelling in opposite directions ; the upgoing wave in 
each layer is the superposition of the reflection of an downgoing wave 

in that layer and the transmission of an upgoing wave from the lower 
layer; and the downgoing wave in each layer is the result of the 

transmission of a downgoing wave from the upper layer and a reflection 



of an upgoing wave in that layer (Fig 15). This concept is similar to 

that for a travelling wave in the acoustic tube , the model for vocal 

tract , in speech processing. Under the modeling constraints ( namely 

normal incidence , one dimensional wave equation and homogeneous & isotropic 
horizontal layers ) these physical facts can be employed to derive 

interrelationships which relate the upgoing and downgoing waveform from 
one layer to the other. These lead to lattice - like structures which 
give physical meaning to the acoustic wave propagation between different 

layers . 

The seismic problem is somewhat different from the speech analysis . 
In seismology , it is not the transmitted wave that is experimentally 
accessible , rather it is the overall reflected wave. On the basis of 
this reflected wave , the layered structure ( i.e the reflection 

coefficients ) are to be extracted using deconvolution. 

However the point to be emphasized is that both the dynamic 
modelling of the earth and the acoustic tube model of the vocal tract 

lead to the " lattice structures ", and hence the analysis is valid for 

both the applications, However it should be noted that the Random model is 

applicable only in seismic signal processing. 

1.6 SCOPE OF THE WORK : 

Each of the two deconvolution methods, viz., Homomorphic and 

Predictive , has particular advantages and limitations when used 

separately. In this thesis , an attempt has been made to combine the two 
methods to make use of the advantages of both. This approach is seen 

to result in significant improvement in the performance. The deconvolution 


method is illustrated by means of synthetic examples for the two models 
viz., Random Model and Dynamic Model. The output of the system , for 
both minimum phase and maximum phase input is synthesised and then it 


is used to 

test 

the 

deconvolution algorithms. First 

the 

Predictive 

Deconvolution 

method 

and 

the Homomorphic Deconvolution 

method 

are used 


separately to recover the convolutional components. The results thus 
obtained are then compared with those obtained after combining both the 

methods. For the Dynamic model we then consider the recovery of the 

system parameters using System Identification Method and Dynamic Predictive 
Deconvolution method. The system identification method involves solution of 
Yule-Walker equations. However this method requires accurate a-priori 
knowledge of the order of the ARMA model of the system. The effect 

of additve white guassian noise on the deconvolution process is also 
studied. The Homomorphic signal processing is found to be more suitable 

for the recovery of the input signal to the system in the noisy as well 
as noiseless cases. 


1.7 ORGANISATION OF THE THESIS 


The 

Thesis 

is divided 

into 

six 

chapters. 

In 

Chapter - 2, 

the 

computation 

of the 

synthetic 

trace 

, for 

the two 

models , 

Random 

and 

Dynamic , 

is discussed. The 

two 

models 

, as applied 

to 

speech 

and 


seismic signal processing is discussed. The modeling of a layered medium 
as a lattice filter is next discussed. 

In Chapter-3, the two deconvolution methods viz., predictive and 
homomorphic are discussed. The deconvolution of the synthetic traces 
using the two methods is illustrated. The problems associated with each 


method are studied. The predictive deconvolution involves the solution of 
a set of normal equations. Levinson's algorithm is used to solve the 
normal equations. A strategy to determine the exponential weighting of 
the synthetic trace, so as to make the system impulse response a minimum 
phase sequence is proposed. 

In Chapter-4, it is shown that there is a significant improvement in 
the results when the two methods — predictive and homomorphic — are 
combined. The all -pass filter for a mixed phase signal which is a major 
requirement in the predictive deconvolution method is determined using 
Homomorphic filtering. Both homomorphic filtering and the Dynamic 
Predictive Deconvolution method together are used to deconvolve the 
synthetic trace for the dynamic model. 

The influence of additive white noise on homomorphic filtering is 
studied in Chapter-5. The thesis is concluded in Chapter - 6, which gives 
a summary of the results obtained . 


The properties of the complex cepstrum are mentioned in Appendix. 



FfGj 1-3 



(i-+) C- L+ +) ’ 





t7 CO U - •!“ 


T1 

S' 


-n 

«1 


m 

H 

m 

n 


Ol 

o 

c 

P 

O 


-4 m 

o 

<o 


Z 

s 

m 

t> 

TJ 

m 

» 


z 

p 

rp 

< 

m 

P 

o 

m 

p 

> 

H 

Z 

$ 


l 

> 

< 

m 



P 

> 

o 

m 

o> 


H 

*» 

a* 

o 

m 

05 


-4 

JU 

> 

n 


m 

0) 


f w ^ 

4o *»' * 
* K ^ 


X X 



3> > 


8 

2 

2 

i 


m 

■o 

-i 

x 


A 

0 

2 

2 

5 

z 

fc» 

m 

■D 

H 

X 


3 ? 
2 z 


I 

J> 

< 

m 

> 


o 

2 

2 

0 

2 


M 

X 

"C> 

O 

Z 

-4 



CHAPTER - 2 


SIMULATION 

2.1 COMPUTATION OF SYNTHETIC TRACE 

A synthetic trace is the output of a known system for a known 

input. Computation of a synthetic trace requires the modeling of the 

system. After assuming some model for the system, the synthetic trace is 
computed by convolving a known input with the impulse response of the 

modeled system. This synthetic trace is then subjected to the 

deconvolution process. The deconvolution algorithms try to recover the 
input signal and the parameters of the modeled system from the synthetic 
trace. The relative performance of the deconvolution algorithms is 
studied by comparing the recover'd signals/parameters with those used 

for the computation of the synthetic trace. These simulation results are 
used to illustrate the improvement in the results of the deconvolution 

process when both predictive and homomorphic deconvolution methods are 
used together. 

Two models for the system are considered and the synthetic 
traces are computed for these models. The models that are considerd are 
equivalent to those considerd for the vocal tract in speech analysis and 
the earth in seismic signal processing. But the parameters of the model 
that are considerd for the computation of the synthetic trace may not 

actually represent either the earth or the vocal tract. 

The modeling of the system and the computation of the synthetic 
traces are discdussed in this chapter. 


2.2 RANDOM MODEL 


The synthetic trace t(n) is the convolution of three signals , and 
is given by , 

t <n) = s <n)* r <n)* p (n) (2.1) 

where , s (n>, r (n), p <n) are all real causal and stable. 

The component s(n) is either a minimum phase or a mixed phase 
signal. The sequence r (n) is a random sequence C eq. 1.3 3. The third 

component p (n) is a periodic sequence given by, 

p <n)- « k tS(n-kN)3 , <2.2> 

k 

where , la I < 1 and N is a positive integer and k « 0, 1, 2, 

The above model is applicable in seismic signal processing , where , 
the seismogram is modeled as a convolution of the seismic wavelet (which 
is either minimum phase or mixed phase), the reflection coefficients of 

the layers of the earth and a reverberation wavelet ( corresponding to 
s(n), r (n) and p (n), in the above model ). The reverberation wavelet is 

present when the source is buried under a strong reflector e.g. water in 

marine seismology. The reverberation effect will occur as a result of 
the trapping of most of the source energy in the water layer. The 
reflectin coefficients of the layers of the earth can be assumed to be 
random. 

2.2.1COMPUTATION OF THE SYNTHETIC TRACE FOR THE 

RANDOM MODEL 

The computation of the synthetic trace is carried out in two 
stages. First the component signals s(n ), r (n ) and p(n) are generated. 



Then all the three components are convolved to get t <n ). 


The signal s(n) cam be either minimum phase or 
The minimum phase signal is represented by s 1Bin (n) and it 
have the following z-transform [173 




(1 - rA" 1 ) (1 - r.e A” 1 ) 


maximum phase . 
is assumed to 

(2.3) 


The mixed phase signal 
is taken to be equal to, 


where , rs 0.9, © = ^ . 
is represented by s lftlx (n) and 


S mto (z) 


( 1 - r/z' 1 ) ( 1 - r.e" je z” i ) 


its 


z— transform 


(2.4) 


where , r = 0.9 , 0= 2 . 

D 

S m# x<z>- <— 0.5— j 0.5— z —i ) (-0.5+j 0.5-2 (0.83-z -1 ) (0.834z“ 1 )(0.91-z“ 1 ) 

- 0 .313 + 0.283 z -1 - 0.516 z“ 2 - i.i z" 3 + 0.09 z“ 4 + z~ s (2.5) 

The signals s min <n) and s miJC (n) , computed using the above z— 

transforms are shown in Fig 2.1 and Fig 2 2 respectively. 

Figs. 2.2 Cal and 2.2Cb3 show the minimum phase equivalent signal and 
the corresponding all-pass filter response , the details about which are 
given later in Chapter-3. The point , which is to be noted at this stage, 
is that the convolution of these two signals will give the mixed phase 
signal s miJC (n) . 



TABLE - 1 



n 

r(n) 

n 


r(n) 


0 

Q .4300 

200 


-0.1765 


14 

-0.0730 

203 


-0.1089 


53 

-0.2930 

241 


-0.1559 


80 

-0.1427 

247 


0.0150 


105 

-0.1465 

275 


-0.2028 


110 

0.2366 

317 


0 .0300 


115 

0.2100 

332 


0.1277 


162 

0.1277 

354 


-0.0638 


167 

0.1500 

377 


-0.1540 


187 

-0.0500 

381 


0.1145 


194 

0.1070 

431 


-0.0526 . 

The 

above 

random sequence 

is a modified 

version of the one given 

in C1JJ. 

The sequence r(n) is 

shown in Fig 

. Z.3. 

The autocorrelation 

function 

*rr<n) 

of r (n) for n £ 

0 is shown in 

the 

Fig 2.5. It can be 

seen that 

/ 







1 frrln) 1 < 

0.1 ^rr(O) for n #0. 


So rln) can be 

considered to be 

approximately 

white , 



The sequence pin) is taken to be the following sequence : 

p(n) = X<-0.387 • 5<n ~ kN) (2.6) 

k=o 

The maximum value of "k" is taken to be 13 because the value of 
(0.387) k becomes comparatively- neglegible for k>13. The Figs. 2.4 and 2.9 



show the p(n) for the two values of N : N=22 and N=40. 


The synthetic trace is computed by convolving the three components 
s<n), r(n), pin >. Figs. 2.7, 2.8, 2.11 show the synthetic trace computed for 
three different combinations of the three component signals. 

Fig. 2.7 ti<n)- Convolution of s m!n (n), r(n) and pin) for N= 22. 

Fig. 2.8 t a (n)- Convolution of s mix (n), rin) and pin) for N = 22. 

Fig. 2.11 : t 3 (n)~ Convolution of s min (n), rin) and pin) for N =* 40. 

2.3 DYNAMIC MODEL : 

A lossless layered medium which is described by the wave-equations 
and boundary conditions are considered for this Dynamic model. Specific 
applications of such layered media are : (1). horizontally stratified 

nonabsorptive earth with vertically travelling i.e. normal incident, plane 
compressional wave, (2). interconnection of lossless transmission lines, (3). the 
vocal tract modeling in speech analysis, etc. 

A layered system consisting of K layers is shown in fig. 2.15, The 
thickness of each layer is considered to be equal. The input signal s(n) 
is applied at the interface #0. Let C; be the normal incidence reflection 
coefficient associated with the j-th interface. Reflection coefficient c$ 
depends on the impedance of the materials at the ’ j ' th interface and 
I cj I < 1. The output y(n), which is observed at the interface #0 contains 
information about the reflection coefficients. Since the characteristics of 

the layers are unknown , the main objective is to reoover the reflection 

coefficients at each interface , which in turn are related to the 

impedances of the adjoining layers of the interface , by the relation , 


where Z; is the impedance of "j" th layer. 

The model assumes normal incidence ; horizontal layers , where each 
layer is lossless, homogeneous and Isotropic ; small strains ; and an one 
dimensional wave equation. Each layer is assumed to have the same 
one-way travel time for a wave propagation. This common one-way travel 


time in each layer is 

taken 

to be one 

half 

a unit time . 

The layers 

having different travel 

times 

are built 

up 

by inserting 

hypothetical 


layers , whose reflection coefficients are equal to zero. 

The upgoing wave in each layer is the superposition of the 

reflection of a downgoing wave in that layer and the transmission of an 
upgoing wave from the lower layer . The downgoing wave in each layer is 

the result of a transmission of a downgoing wave from the upper layer 
and a reflection of an upgoing wave in that layer [ Fig. 2.16 1 Let d/n) 

and u/n) be the downgoing and upgoing waves respectively in the j-th 

layer. The waves at the interface "j" can be related to those in the 

(j+i) th layer 1173 as follows : 

(i-Cj) . dj+ifn) = dj(n— | ) — Cj . Uj(n+|) ( 2 . 8 ) 

(i — Cj) . u j+ 1 (n) m Uj(n+|) - Cj . dj(n— | ) ( 2 . 9 ) 

Using these equations a lattice form can be established which 

implements the upgoing and the downgoing wave in any layer. This lattice 

form is shown in the Fig. 2.17. The above equations can be written in 
matrix form C173 as follows : 


Dj+i(z) 

,1/2 

1 

••s 

0 

1 

1 

N 

1 


Dj(z) 


» _£ 




Uj+ife) 

i-Cj 

-CjZ -1 i 


Uj(z) 



m m . 


. « 


<2.i0> 



Expanding (210) , a relation between the upgoing and the downgoing 


waves in the layer #0 to those in (K+i)th layer can be written as 


■ 


- 


» ■ 

D K+i (z) 


Pk* Qk* 


D 0 (z) 


tt/2 




U«+i(z) 

t k 

Qk p« 


U 0 (z) 

. „ 


. - 


. 


where, 


Pk* GU* 

Q* P* 


K 

n 



( 2 . 12 ) 


and , 




T k = 

K 

n ( i 

- C s ) 

<2.!3) 




$■1 


To 

compute the 

output of 

this layered system for 

a unit impulse 

input , 

the boundary 

condition 

u K+1 (n) 

= 0 is assumed 

to hold. This 

corresponds to the 

case when 

there 

is no upgoing waves from the 

deepest 

layer i.e. the fl<+i) th 

layer 

and implies that 

at the final 


interface , all the downgoing energy is absorbed. 


The layered medium model is also applicable to speech analysis. The 
similarity between the modeling of the layered earth in seimic signal 
processing and the vocal tract in speech analysis is discussed in [43 
which describes the vocal tract lip reflectance by an ARMA filter. 


2.3.1 COMPUTATION OF THE SYNTHETIC TRACE FOR THE 

DYNAMIC MODEL 

The reflection coefficients, c n 's for the dynamic model are taken to 

be such that c n+8a =r(n) , where r(n> , n=0,i, 43i is the sequence given in 

Table— I with oq 4 0.9. Thus two-way travel time in the layer#! is taken 



to be 22 units. The output is considered to be recorded just below the 
interface #0. 

Using (2.12) the polynomials P k and Q k are next computed using the 
above values of c/s . For the boundary conditions U k+ i(z) = 0 and for 
unit impulse input, i,e. D 0 (z)=i , we get from the relation-(2.ii) 

Uoft) = - . D 0 (z) (2.14) 

The sequence Ui(n) is computed from the relation— (2.11). 

Finally , the two synthetic traces t 4 (n) and t £ (n) are computed by 
convolving Ui(n) with the input s min (n) and s mix (n) given in section 2.2.1 

The Figs. 2.12—2.14 show unit impulse response u t (n) and the 
synthetic traces t 4 (n) and t 4 (n>. 


A 

fast 

algorithm 

to compute the 

polynomials 

Pk 

and Qk Is 

suggested 

by 

Choate, C61 , 

which reduces 

the number 

of 

multiplications 


from K 2 to K.( In K f . 



smin 


Fig. 2.1 

Minimum phase signal 










g.2.7 The synthetic trace 

obtained bv convolving 
smin(n), r(n) & p(n) ~ 
period of p(n)=22 


The synthetic trace 
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CHAPTER 3 


PREDICTIVE AND HOMOMORPHIC 
DECONVLUTION METHODS 

The synthetic trace computed in the previous chapter is available 

for deconvolution. It is assumed that none of its convolutional components 
are known ; different deconvolution methods try to recover the 
convolutional components by making suitable assumptions. The recovered 
signals are then compared with the original signals and the efficiency of 

the different deconvolution methods are studied. The inference about the 
relative merits and demerits and the complexities associated with each 
deconvolution method is obtained. In this chapter the procedure is 
carried out for the Random model . 

In this chapter, the predictive deconvolution and the Homomorphic 
deconvolution methods are used separately to recover the convolutional 

components. In chapter— 4, we show that better results can be obtained 

when both the methods are combined together. 

3.1 PREDICTIVE DECONVOLUTION 

Synthetic trace : t( n ) = s( n ) * r( n ) * p( n ) (3.1) 

Assumptions made: 

1. s( n ) is a minimum phase signal. 

2. r( n ) is a white noise sequence. 

Consider the z— transform of the synthetic trace: 


T( z ) as S( z ) . R( z ) . P( z ) 


< 3 . 2 ) 


or, T( z -i ) = S( z~ l ) . R( z -i ) . P( z -1 ) (3.3) 

Multiplying both the equations , 

T( z ).T( z -i ) = S( z ).S( z -1 ) . R< z ).R( z - - 1 ) . P< z ).P( z -1 ) 

Taking the inverse z— transform of both the sides , we get , 

n ) = * s ( n ) * M n ) * * P ( n ) (3.4) 

where , ^ s ( n ) > <M n ) , ^ P ( n ) respectively are the autocorrelation functions 
of s( n ) , r( n ) and p( n ) . 

By assumption -2 , the sequence r( n ) is white , it follows that 

i 

n ) = R.S(O) 

where R is is a constant. 

n ) = R. * s ( n ) * <M n ) (3.5) 

Thus the autocorrelation of the synthetic trace is a scaled version 

of the autocorrelation of the convolution of the two convolutional 
components , s( n ) and p( n ). First let us consider the deconvolution of 
both s( n ) and p( n ) from the synthetic trace , so as to recover the 

sequence , r( n ). That is , it is required to find a filter f( n ) , which 

when convolved with the synthetic trace , produces the sequence r( n ). 

ry( n ) f f( n) = r( n ) * s( n ) * p( n ) * f( n ) = r( n ) (3.6) 

s(n)*p(n)H(n) = 5(n) (3.7) 

The signal, s( n ) X p( n ) is not known. However, a scaled version 

of its autocorrelation , which is the autocorrelation of the synthetic 

trace itself [eqn. 3.5] , can be computed. The scaling does not affect the 

design of the filter f(n). A least— squares filter which minimises the 

error between the desired output and the actual output is designed. The 

expression for the error is given by, 

E - £ t z( n ) - y( n ) f 
i 



where , z(n) = S(n) is the desired output and y( n) = s( n ) * p( n ) # f < n ) 


is the actual output. 

The minimisation of the error E , results in the Wiener Filter, f( n ) , 


which 

can be 

obtained by 

solution of 

the r 

* t <0> 


.... * t <L-l> 





f t (0) 

.... * t <L-2> 


f(0> 


1 

.... 

.... 




f(i> 


□ 

L *t<L-i> * t (L-2) 

o 








f(L-i) 


0 


(3.8) 


where, ^ t (n) is the autocorrelation of the signal s(n)#p(n ) , i.e. of the 
synthetic trace. 

The normal equations are solved using Levinson's recursion algorithm 
which exploits the Toeplitz structure of the autocorrelation matrix. The 

total number of multiplications required to solve for a digital filter with 
L coefficients is proportional to L z for Levinson's recursion method , as 

compared to L 3 for Gauss-elimination method. Levinson's algorithm 

requires computer storage space , proportional to L rather than L 2 as in 
other methods. The filter f(n) obtained from the solution of the above 
normal equations is next convolved with the synthetic trace , t(n) to get 
back the random sequence r(n). 

Let s(n) be a mixed phase signal. Any mixed phase signal can be 

represented as the convolution of its minimum phase equivalent signal and 
a corresponding all-pass filter: s(n) = s^^n) * h 4p (n) , where , s me ,(h) is the 
minimum phase equivalent of s(n), h 4p (n) is the impulse response of the 

all-pass filter. 



So the synthetic trace , t(n) is written as : 

t(n) = r(n) * h ap (n) * s meq (n) * p(n) 

The filter f(n) obtained by solving the eqns. 3.8 is such that, 
f(n> * s maq (n) « p(n) = S(n>. 

So, f(n) * t(n) = r(n) * hip(n) is the recovered r(n) , which contains an 

all-pass component. 

Now the sequence s(n) is to be recovered. s(n) can be either a 

minimum phase or a mixed phase signal. We now show that the s(n) can 

be recovered only if it is minimum phase and the period of the 

sequence p(n) is greater than the length of s(n) . Let N be the period 

of p(n). Let the length of s(n) be less than N. Then the first N samples 

of the signal x(n) =s(n)#p(n) correspond to s(n) itself. Now we design a 

f'(n) which is a prediction operator with prediction distance equal to N. 

The output of the filter f'(n) is given by, 

y’(n) = 5£x(i)*f(n-i) = *(n+N) (3.9) 

i 

where x(n+N) is an estimate of x(n+N). 

The prediction operator is computed by minmising the error E given by, 

E - £ C z'< n ) — y'< n ) f 
i 

where, z'(n)=x(n+N) is the desired output and y’(n) is the actual output. 

The minimisation of E leads to the following matrix equation : 

* t co) * t a> 

*t<0> 


^ t (K-2> 

* t <0> 


f'(0) 


ft<N) 

f'(l) 


*(N+i) 

f'(K-i) 


*<N+K-i> 


(3.10) 



For a prediction operator f'(n) of length K the corresponding prediction 
error operator f"(n) with prediction distance N is : 

f”(n> = C i, 0, 0, A -f'(0), -f'<l), -f'(K-l) 3 <3.ii) 

Convolution of x(n) with f"(n> makes x(n) vanish for the lags greater 
than N . 

Thus if the duration of s(n> is less than the period N of p(n) 

then convolution of x(n) with f"(n) gives s(n). We do not know x(n) 

separately. The prediction error operator is hence convolved with the 
trace t(n). As r(n) is a white noise sequence, the autocorrelation of the 
resultant trace is the autocorrelation ^ s (n) of s(n> itself. Next the AR 
filter coefficients f(n> are computed from eqn. 3.8 with ^ t (n) replaced by 
f s (n). Then s(n) is simply the impulse response of this AR filter. If s(n) 
is a mixed phase signal then s(n) recovered from the above method will 

be its minimum phase equivalent signal. 

The predictive deconvolution fails to recover the signal s(n) if the 
period , N is less than the duration of s(n). If the wavelet is mixed 
phase , then only its minimum phase equivalent signal can be recovered. 

The Predictive deconvolution method requires r(n) to be a perfect random 
sequence , which is not realisable. 

3.1.1 ILLUSTRATION OF PREDICTIVE DECONVOLUTION : 

(a). Consider the synthetic trace , ti(n). 

1 . Find the autocorrelation of ti(n). 

2. Solve the equatios ,3.8 , for f(n) using Levinson's recursion. The 
filter length L is taken to be equal to 600. 

3. Convolve f(n) with t*(n) to recover r(n). The recovered r(n) is shown 



in Fig. 3.1 2 along with its actual values [ from Fig 23 ] . The deviation 
of the recovered r(n) from its actual values is due to the fact that 

r(n) is not a perfect white noise sequence. Hence in the above 

procedure, the minimum phase equivalent of r(n) also gets deconvolved 
along with s(n) and p(n). 

4. Assume that the duration of s(n) is approximately 40 units Esee 

Fig2.il. Use the method explained above to recover s(n). The recovered 
s(n) is shown in Fig. 3.1.1. The degradation of s(n) is due to the 
period of p(n) being less than 40 ( equal to 22 for trace ti<n) ). 

(b>. Consider the synthetic trace , t a (n) 

The procedure given in (a) is repeated. The recovered s(n) and 
r(n) are shown in Figs. 3.1.3 and 3.1.4 , respectively. Both the recovered 

signals deviate considerably from the actual signals, because s(n) is a 

mixed phase signal s* iK (n) in the trace t 2 (n). 

<c>. Consider the synthetic trace , t 3 (n). 

The autocorrelation of t 3 (n) is shown in Fig.3.1.5(a). The same 

procedure as in (a) is repeated for ta(n). The recovered s(n) is shown 


in Fig. 3.15 

This recovered s(n) is 

a 

very 

good 

replica 

of the actual 

s(n), because 

in this case N=40. 






Thus 

the problems which are 

not 

tackled 

by 

the Predictive 

Deconvolution 

are the determination 

of 

the 

all-pass 

filter 

when s(n) is 

mixed phase 

, and the restoration 

of 

the 

minimum 

phase 

equivalent of 


r(n), which gets deconvolved when r(n) is not a perfectly random 


sequence. 



3.2 HOMOMORPHIC DECONVOLUTION 


Homomorphic decconvolution does not assume any of the convolved 
components to be minimum phase or perfectly random. However the quality 
of the results depend upon the degree of separatability of the component 
signals in the " quefrency " domain. Two convolutional components are well 
separated in the quefrency domain if the spectrum of one component is 
slowly varying or smooth and the other is rapidly varying. The efficacy 
of this method depends mainly on factors such as the duration of s(n), 


the 

time 

difference between 

the 

first 

two 

nonzero 

samples 

of 

the 

sequence 

r(n) , the 

period of 

P<n) , 

the 

choice 

of the 

window width 

etc. 

The 

determination 

of the amount 

of 

the 

exponential 

weighting 

of 

the 


synthetic trace in order to make the sequence r(n) , minimum phase also 
plays an important role in the success of the Homomorphic deconvolution 


method. A trial 

& error 

method for the 

determination 

of 

the 

amount 

of 

the 

exponential 

weighting 

is proposed 

in 

this thesis. 

A 

brief 

review 

of 

the 

Homomorphic 

Signal 

Processing 

is 

given in 

the 

Appendix. 

The 


deoonvolution method is explained below. 

Consider the synthetic trace , t(n)= s<n) * r(n) * p(n) 

The spectrum of s(n) is smoothly varying and that of r(n) and p(n) 

are rapidly varying. So the complex cepstrum of s(n) tends to be 
concentrated at the low quefrencies. If the signal s(n) is a mixed phase 

signal , its complex cepstrum is nonzero for n<0 . The sequence p(n) is a 
minimum phase periodic sequence with period N. Hence its complex 
cepstrum is zero for n< N and is a periodic sequence sequence with the 

same period N. The random sequence r(n) in practice is a mixed phase 



signal. So its complex cepstrum occupies the entire quefrency domain. The 
complex cepstrum of r<n) can however, be made zero for n<0 by a 
suitable exponential weighting of the synthetic trace , such that r(n) now 


becomes a 

minimum 

phase 

sequence. Let the 

synthetic trace 

be 

exponentially 

weighted 

by an 

amount, "a" , [ a<i ]. 

The z -transform 

of 


a n .t(n) is 



T(az —i ) = 

S(az -1 ) . R(az -i ) . P(az _1 ) 


(3.12) 

Thus if the 

synthetic trace 

is exponentially weighted , then 

each 

of its 

convolutional 

components also 

gets exponentially weighted 

by the 

same 

amount. Since 

a<l , all the 

zeros of the R(az — *) can be 

made 

to lie 


within the unit circle by a proper choice of "a" , and thus r(n) can be 

made minimum phase. Then the complex cepstrum of this minimum phase r(n) 
is zero for n< M+i , where , M is the the number of samples between 

the first two nonzero samples of r(n). The complex cepstrum t(n) , of the 
exponentially weighted synthetic trace is next computed. Now , only the 

complex cepstrum of s(n) is present near n=0. So s(n) can be obtained 

A 

low— pass "Uttering" of the complex cepstrum , i.e. by multiplying t(n) by 
a window , w(n) , given by : 

w(n) = 1 , n< L ; 

= 0 , elsewhere. 

where , L = min(N, H+l). If s(n) is negligibly small for n> L , then s(n) 
can be recovered successfully from s(n>, no matter whether it is minimum 
phase or mixed phase. The signal s(n) thus recovered is to be 

exponentially deweighted at the final stage of the reconstruction. 

The inverse filter f(n) , for the deconvolution of s(n) is obtained in 
a similar manner by using a window , ww(n) = — w(n) , where w(n) is defined 



above. The filter f(n) is exponentially deweighted and then convolved with 
the synthetic trace to deconvolve s(n). The signal s(n) can't be 

deconvolved using high— pass liftering because exponential deweighting will 
amplify the noise due to aliasing in the quefrency domain at the end of 
the deconvolved trace. This aliasing in the cepstral domain is always 
present because the complex cepstrum of a finite signal is of infinite 

duration. The aliasing in the quefrency domain can be reduced by taking 
the FFT and the IFFT length equal to about four times the length of 
the synthetic trace . 

It is not possible to deconvolve p(n) from the trace, because the 

complex cepstra of r(n> and p<n) overlap. However, the sequence p<n) can 
be suppressed in the particular case in which the first few cepstral 

components of p(n) do not overlap with the cepstrum of p(n) and an 
a-priori knowledge of the period of p(n) is available. {For example, 
consider the complex cepstrum of the appropriately exponentially weighted 
synthetic trace t t (n). The first two nonzero cepstral components of r(n) 
are present at n=14 and n=53 whereas those of p(n) are at n=22 and 

A 

n=44 }. Then the two cepstral values of t(n) at quefrency, ns= N and n= 2N 
are made zero, This eliminates the first two samples of p(n) at n= N and 
n= 2N & reduces the remaining values of p(n) by i/3 of their original 
value. Next, without any low-pass liftering , the synthetic trace is 
reconstructed after exponential deweighting. The exponential deweighting is 
done for only half the trace in order to reduce the amplification of 
noise due to aliasing. The contribution of p(n) is insignificant in this 
reconstructed trace . The inverse filter f(n) as computed above is then 
convolved with this reconstructed trace to recover the sequence , r(n). 



The 

determination of 

the 

amount of 

exponential 

weighting 

and 

the 

width of 

the window are 

the 

two main 

problems associated 

with 

the 

homomorphic 

deconvolution. 

So 

far , no 

quantitative 

method 

for 

the 


determination of the amount of the exponential weighting is proposed C241. 
In this thesis , a trial & error method for the determination of the 
required exponential weighting is proposed, which gives satisfactory 
results. This method is explained in the next section. The window-width 
is decided through observation. 


3.2.1 DETERMINATION OF EXPONENTIAL 'WEIGHTING FACTOR; 

Let "a" , where a<l , be the amount of the exponential weighting 
which is sufficient to make the sequence r(n) a minimum phase sequence. 


Then both the real and the 

complex 

cepstra 

of 

the sequence, 

a*. 

r(n) 

will 

be zero around 

zero and 

negative 

quefrencies 

and 

will 

be 

equal 

at 

higher positive 

quefrencies . 

This 

property 

is 

made 

use 

of 

in 

this 


method. Here, the high quefrency parts of the real cepstrum and the 
complex cepstrum of the synthetic trace , are compared for different 
amounts of exponential weighting , and the amount for which both the 
parts are same is chosen as the required exponential weighting. The 
values of "a" are chosen in a decreasing order starting from a=i. The 
oepstrum of the s(n) is concentrated only around zero quefrencies & the 
cepstrum at higher quefrencies is mainly due to r(n) and p(n). Hence this 
method works quite satisfactorily. 



3.2.2 ILLUSTRATION OF HOMOMORPHIC DECONVOLUTION 
Consider the synthetic trace , t t (n) . 

1. Use the trial & error method to determine the exponential 

weighting required. The complex and real capstra of a n .t(n) (a=0.992, 0.989) 
are shown in Fig 32.1(a) - 322(b). It can be seen that the high 

quefrency values of the complex cepstrum and the real cepstrum are the 
same for a=0.989 and are not same for a =0.992. So the required 
exponential weighting is taken to be 0.989. The complex cepstrum and the 
real cepstrum computed for this exponential weighting are used for 

deconvolution. 

2. In Fig.322(b), it can be seen that a strong peak appears at 

n=14, as the first cepstral value of r(n) . So the window w(n) is fixed 

through this observation , to be of length L=i4. 

3. The s(n) recovered from the real cepstrum is shown in Fig. 3 23. 
This recovered s(n) from the complex cepstrum of the synthetic trace will 
also be the same because the s(n) in the synthetic trace t A (n) is 
minimum phase. 

4. The inverse filter f(n) of s(n) is computed by multiplying the 
cepstrum by a window ww(n)=— w(n). 

5. Assuming that the period of p(n) is known a-priori to be equal 

to 22, the cepstral values for n=22 and n=44 are made zero. This will 
suppress remaining values of p(n) by about i/3 . The synthetic trace is 
then computed back. It is then convolved with f(n) to get back r(n). 

The recovered r(n) is shown in Fig 32.4. 



Consider the synthetic trace , t 2 (n). 

The same procedure is repeated for this trace also. The complex 

and the real cepstra of a n .t 2 (n) (a =0.989) are shown in Figs. 325(a) and 
325(b). The minimum phase component of s(n) is recovered (Fig 3 2. 6) from 
the complex cepstrum by using a window, w(n) given by, 
w(n) = i , n< L ; 

= 0 , elsewhere. 

The maximum phase component of s(n> obtained from the negative 
quefrencies is shown in Fig. 32.7. Convolution of these two components 
gives the mixed phase s(n) ( Fig. 32 B ). The s(n) recovered from the real 

cepstrum is shown in Fig. 3.2.9. This is the minimum phase equivalent of 
the Smjxfn). 

Although the homomorphic deconvolution does not make use of any 

of the basic assumptions regarding the phase of s(n) and the randomness 
of r(n) , nevertheless it is a deterministic processing method. Hence , 
there is more uncertainly about the quality of the results obtained 
through this deterministic procedure. Choice of the width of the window 
and the amount of exponential weighting should be quite accurate. Slight 
deviation from the exact value of these two parameters yields poor 


results . 
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CHAPTER 4 


DECONVOLUTION BY COMBINED PREDICTIVE 
AND HOMOMORPHIC FILTERING 

As pointed out in Chapter-3, the problems associated with the 
Predictive deconvolution method are : (1). determination of the allpass 

filter when s(n) is minimum phase , (2). poor results when the sequence 

r(n) is not a perect random sequence. Also the predictive deconvolution 
method requires an estimate of the duration of s(n). As discussed in 

section 3.2 , the Homomorphic deconvolution method in general can not 
recover r(n). In this section , it is shown that if both Predictive and 

Homomorphic deconvolution methods are combined , retaining the assumptions 
made by the former , a significant improvement in the quality of the 
results can be obtained. Through this method , the problems associated 

with Predictive and Homomorphic methods when used independently, are 
solved. In this method , the all-pass component of s(n) and the minimum 
phase equivalent of r(n) are computed using homomorphic filtering. 

4.1 DETERMINATION OF THE ALL-PASS FILTER 

A stable and causal minimum phase sequence has a property that 
all the zeros of its z-transform lie within the unit circle. But a 
stable and causal mixed phase signal will have some of its zeros 
outside ■ the unit circle. This mixed phase sequence can be represented 
by a cascade of a minimum phase system and an all-pass system having 



unit magnitude frequency response. A first order all -pass filter is 


represented by , 

H 4p (z) = — — , 0<a<i , a real. (4.1) 

(l-az -1 ) 

Consider a mixed phase system X(z) , with one zero outside the 
unit oircle, at z<* ^ , ( la I < 1 , "a" real ), and remaining zeros inside the 
unit circle, i.e. 

X(z) - Xl(z) . <z -i - a) (4.2) 

where , Xl(z) is minimum phase. 

Equation 4.2 can be written as , 

X(z) = Xl(z) . (z _1 - a) . (i ~ a - z ~. 1 - ) 

(1-a.z A ) 


= [ Xi(z) . (i-az -1 ) ] . t l ~A 3 

(i-az ) 

- X n e*(z) . H 4p (z) (4.3) 

Xj.(z) is called the minimum phase component of X(z) and X m6l ,(z) is called 
the minimum phase equivalent of X(z). I H ap (z) I = i, for all frequencies. 


In eqn. 4.2, x(n) can be considered as the convolution of its 
minimum phase component x 4 (n) and its maximum phase component x 2 (n) , 


where, 

X 2 (z)= z 1 — a. 

In 

general, if 

X(z) has P zeros outside 

the unit 

circle, 

then it can 

be 

expressed 

interms of its minimum and 

maximum 


phase components as-' 


X(z)=Xi(z).X 2 (z) 

p 

where, X 2 (z)=l+^ a-, z~' 

i=i 

As in (4.3) X(z) can also be expressed as: 

X(z)=X me «,(z) .H ap (z) 


where H ap (z) is of the form, 
H ap (z)= 


+ # a z 


A + + /3pZ ** 


£p + £p—i,z 1 + 


+ 0 x z p 1 + 1 


(4.4) 

(4.5) 


(4.6) 



Hence it follows that all the zeros of H ap (z) must be the same as those 

of X 2 (z), which implies that fossa; , 1=1,2, P. Thus we have seen that 

only the knowledge of x 2 (n) is sufficient to determine the corresponding 
all-pass filter. 


Considering the first problem associated with the predictive 
deconvolution , viz., the determination of the all-pass filter when s(n) is 
mixed phase , it is sufficient if we know the maximum phase component of 
s(n) . Homomorphic filtering can be used to get this maximum phase 

component of s(n) . The random sequence r(n) is made minimum phase after 
exponential weighting of the synthetic trace and then the complex 

A 

cepstrum t(n> is computed. It is assumed that the exponential weighting 
does not shift any of the zeros of sCnl lying outside the unit circle to 
inside the unit circle . The maximum phase component of s(n) whose 
complex cepstrum is present only at the negetive quefrencies can be 

A 

recovered by considering the samples of t(n> for n<0. This recovered 
signal is just the sequence, { 1, «i, a Sl a P }. Thus , the required all- 

pass filter transfer function is determined. The impulse response of this 
all-pass filter is next computed. The predictive deconvolution requires 
the deconvolution of this all-pass filter from the recovered r(n). The 
impulse response of the inverse of this all-pass filter is just the time 
reversed impulse response of the all-pass filter itself. This is 


illustrated 


for a first order all-pass filter 


H ap (z) 


<z~ x - a) 
(l-az _1 ) 



Me have, 



<4 .7) 


[*«*<*>]-* - fcS, 

- Hapd -1 ). 

So the the time reversed impulse response of the all-pass filter 
is convolved with the the recovered r(n) in Predictive deconvolution 

method to remove the all-pass filter component from it. 

\ 

Treitel S. and E. A .Robinson C22D have suggested a different method 
for the determination of all-pass component . In this method, the mixed 
phase signal s wix (n) and its minimum phase equivalent, s meq (n) are recovered 
from the complex and the real cepstra of a n .t(n) respectively. Next, the 
Smaq<n) is deconvolved from s m!st (n) to get the all-pass component. This 

method requires the accurate choice of the window width. So he has 

concluded that the determination of the window width in Homomorphic 

method is equivalent to the determination of the all-pass filter in the 

Predictive method. On the other hand, the method of determination of the 

all-pass component from the maximum phase component of s TO te(n) does not 
use windowing. This method assumes that even after exponential weighting 
all the zeros of the maximum phase component of s mlx <n) are still outside 

the unit circle. This assumption can be made because the zeros of r(n) 

are much close to the unit as compared to those of s(n). However, we 

shall illustrate in section— 4.4 that, even if few zeros of maximum phase 
component of s m ; x (n) lying close to the unit circle are shifted inside the 
unit circle, the all-pass component so determined does not differ much 

from the true all-pass component. 



4.2 


FINDING OUT THE PERIOD OF p(n) : 


If the Homomorphic Deconvolution method makes use of the 
assumption about r(n) to be nearly random, the period of p(n) can be 

determined. Let r(n) be nearly a random sequence. Hence its autocorrlation 
function has very small values for n?tO. Further, the autocorrelation of 

the periodic sequence p(n), is also a periodic sequence with the same 
period as p(n). Hence the autocorrelation of p(n) has relatively large 
amplitudes for the first few periods. Also , the autocorrelation of the 

convolution of the signals is the convolution of their autocorrelations. 
The real cepstrum of the autocorrelation of a signal is twice the real 
cepstrum of the signal. So the real cepstrum of the synthetic trace is 

computed. This real cepstrum corresponds to a scaled ( scaling factor = 
i/2 ) version of the cepstrum of the autocorrelation of the synthetic 
trace. The large amplitudes at the first few higher quefrencies are 
mainly due to the p(n) . As the autocorrelation of p(n) . is periodic, the 
value of the quefrency , n , at which a large spike appears for the first 
time gives the period of p(n). If the first two cepstral values at n=N , 
& 2N are suppressed , then the convolutional component p(n) is suppressed 

at n=N & 2N and its amplitudes at other n will be i/3 of their original 
values. 

4.3 RESTORING THE MINIMUM PHASE COMPONENT OF r(n): 

If the sequence r(n) is not a perfect random sequence then in 
predictive deconvolution, the prediction error filter computed will be the 

inverse filter of the convolution of minimum phase equivalent of r(n), 
the p(n) and the minimum phase equivalent of s(n). Hence the minimum 



phase equivalent of the r(n) will also get deconvolved In case of 

perfectly random r(n) , its minimum phase equivalent will be a spike at 
nssO and hence its deconvolution does not make any difference If r(n) is 
not a perfectly random sequence , then its minimum phase equivalent will 
be nonzero for n>0 and hence its deconvolution makes a considerable 
difference in the recovered r(n). If this minimum phase equivalent is 
restored back after predictive deconvolution , then the recovered r(n) will 
be correct. This restoration is done as follows : 

The real cepstrum of the synthetic trace , t(n) is computed The 
p(n) is suppressed as explained in section 4.2. The synthetic trace is 
reconstructed . This reconstructed trace is the convolution of the minimum 
phase equivalents of s(n) and r(n) . The minimum phase equivalent of s(n) 

can be recovered from the low quefrency part of the real cepstrum of 
the exponentially weighted synthetic trace. This is deconvolved from the 

above reconstructed trace to get the minimum phase equivalent of r(n) . 

4.4 ILLUSTRATION 

This illustration is to show the improvements in the results obtained 
by the Predictive deconvolution , after making use of the Homomorphic 
Signal Processing . 

,'V.' ••. Consider the synthetic trace , t 2 (i). 

1. The recovery of r(n) using Predictive deconvolution is 

explained in section 3.1 and the recovered r(n) is 
shown in Fig 3.1.4 

2. Recovery of minimum phase eqi valent of r(n) : 


i. The real cepstrum of t 2 (n) is computed (Fig 4.4.1) 



0.989 shifts this zero to inside the unit circle. The mixed phase signal 
sW n >> having this zero in addition to those of s rt , ix (n) is shown in Fig 
4.4.7. The , ;l-pass component corresponding to the zero at z =— i /0992 is 
shown in Fig.4.45. Let s' ttix (n)=s' me ^n).h' Jip (n) as in eqn(4.3). The h' 4p (n) is 
shown in Fig.4.4.9. Then let us try to recover the all-pass component 
h' ap (n> from (0.9B9) r '.s' l#i3C Cn). The complex cepstrum of (0.989> n .s , mix (n) i s shown 
in Fig. 4.4.10. The all-pass component computed from the negative 
quefrencies is shown in Fig. 4.4.11. The all-pass component corresponding 
to the zero at z= -1/0.992 is not present in Fig. 4.4.11. However the 
recovered (Fig.4.4.ii) and the true (Fig.4.4.9) all-pass components are almost 
the same. The reason for this is that the all-pass component 
corresponding to a zero close to the unit circle has its samples for 
n>0 much smaller than the sample at n=0 (Fig .4 .4 .8). 

4.5 DECONVOLUTION OF SYNTHETIC TRACE FOR 

DYNAMIC MODEL 

The synthetic trace for the Dynamic Model i s ; 

t(n) = s(n) * u t (n) 

Both Homomorphic and the Predictive filtering are used to recover 


the reflection 

coefficients 

of 

the layers 

from 

t(n) 

The 

Homomorphic 

Signal 

Processing is 

used 

to 

deconvolve/recover 

the 

input 

s(n) to the 

layered 

system. 

After 

the 

deconvolution of 

the 

s(n) from 

the synthetic 

trace , 

we get 

the 

the 

unit 

impulse response 

Ui(n) , of 

the layered 


medium . A method called Dynamic Predictive Deconvolution Method, proposed 
by Robinson £181 , is made use of to recover the reflectin coefficients of 
each layer from u^n). The procedure for the recovery of s(n) is same 



as the one explained in section 3.2. So only the Dynamic Predictive 

Deconvolution Method is explained in the following section: 


The method is based on the assumption that there is no internal 
loss of energy by absorption within the layers. So the energy from a 
downgoing unit spike at interface #0 as input is divided between the wave 
transmitted by the layered system into the last layer and the wave 
reflected by the layered system above the interface #0. This reflected 
wave is u 0 (n). The reflection (R) and the transmission (T) responses of the 
layered medium are defined as the responses of the medium to an 
impulse incident on the surface as shown below : 


1 R 

i t 

I LAYERED MEDIUM 


i 

T 


V 

t 

LAYERED MEDIUM 
t I 

1 R' 


The relation between arbitrary incoming and outgoing waves , shown 
below , is expressed as 


d 0 (n) 

I 


u 0 (n) 

t 




1 

d K+1 (n) 


t 

u K+i(n) 


From (2.11), 


d k+i(z) 

g K/a 

Pk R 

Qk R 

D 0 (z) 

Uk+i(z> 

T k 

Qk 

Pk 

U 0 (z) 

> ■ 


where , T = XT (i-c s ) , C; : reflection coefficient of ith interface 
s=i 


From the above figures , 

1 

R 


r K/S 

m 

-rK 


Pk 

Qk 


GW 

Pk 


R 


T 

0 


(4.8) 



(4.9) 


0 

T* 



Qk R 



R' 


Pk 


i 


The above two equations can be combined as , 


i 

R 


0 

T' 


_K/a 


P K 
Qk 


Qk 

Pk 


T 

0 


R' 

i 


(4.10) 


Solving , for reflection and transmission responses , we get , 

**> - -?“sr • T(z) - 


- T K . 2 - K/2 


R'(z) = 


Pk<Z) 


Pk(z> 


t K — K/2 

T'Cz) = — L - z 


P K (z) 


Since on physical grounds , the transmission response T(z) must be the 
z-transform of a stable and causal sequence , it follows necessarily that 
P K (z) must be a minimum phase polynomial. 


By the conservation of energy , the input energy is equal to the 
sum of reflection energy and the transmission energy. The input energy 
is unity for unit impulse input. The reflection energy is : R(z)R(z — i ). 


Now, 


1- R(z)R(z -1 ) = 


P^.PkU -1 ) - Q^.QkCz" 1 ) 
PkCzJ.PkCz-*) 


It can be shown that , 

- PkCzI.PkCz^ 1 ) - Q k (z) .Q k (z — 1 ) * f'.T/ 

nr K 

where , T x = TT d+c s ) 

s=l 

So 1 - R(z).R(z —i ) = . T(z).T(z -i ) 

T k 


(4.11) 


(4.12) 


r(n) = c 0 . S(n) + (i- c 0 ).Ui(n) 


(4.13) 



From eon 4 12 , the autocorrelation of the transmission response can be 
computed Since P K (z) is a minimum phase polynomial , we can compute the 

prediction-error operator that contracts the transmitted wave to a spike 

By convolving this prediction error operator with the reflection response. 
r(n) , we get q^n) The coefficients of the polynomial Qk(z) are 
approximately equal to the reflection coefficients. But if the exact value 
of the reflection coeffiocients are required then the following inverse 
recursion formulas will have to be solved. 

q n _ A (n) — (i— On 2 ) - *.[ q n (n> - o n .Pn(n) ] <4. 14) 

P n -i(n) = (1- c n 2 ) — i [ Pn(n) - c n q^n) ] (4.15) 

4.5.1 ILLUSTRATION OF DYNAMIC PREDICTIVE DECON VLUTION 

METHOD 

The synthetic trace is : t(n) = s(n) * u t (n). 

EM. Consider the synthetic trace , t 4 (n). 

1. The minimum phase s(n) is recovered using Homomorphic Signal 

Processing as explained in section 3.2 . The exponential weighting 
required to make r(n) minimum phase is found to be equal to 0.991. 

2. The s(n) is deconvolved from t 4 (n! to get Ui(n) The Uj,(n) obtained 

after the deconvolution of s(n) is shown in Fig 4.5.1. 

3. The reflection response , r(n) is computed from eqn. 4.13 with Co=0.9 

4. The autocorrelation 0(n) of r(n) is computed. 

5. The autocorrelation of the unit impulse input is also an unit impulse. 

So the difference between the autocorrelations of the source signal and 
of the r(n) is an autocorrelation function # , given by , 

*( 0 ) = 1 - 0 ( 0 ) 

0(n) ss -0(n) , for all "n" not equal to 0 . 


& 


Next the prediction-error operator , d 0 =i, d 1( d 2 , .. „,d K <K=500) is to 
be computed by solving the normal equations , 



7. Compute q(n), by convolving the prediction error operator with r(n). 

This is shown in Fig 4.5 2 . 

8. The exact reflection coefficients are obtained by , the inverse 

recursion , given by eqns. 4.14 and 4 15. The reflection coefficients so 

recovered are shown in Fig 4.5.3. 

C II 1 . Consider the synthetic trace t 5 (n). 

In this case , the s(n) is mixed phase . To recover s(n) , then the 
complex cepstrum of the exponentially weighted (a=0.99l) t 5 (n) is to be 
computed. After deconvolving s(n) , the next procedure is same as that 

for t*(n) , given above. 

4.5.2 COMPUTATION OF THE POLY MOM I ALL , Q K (Z) , SOLVING 
MODIFIED YULE-WALKER EQUATIONS 

The Dynamic Predictive Deconvolution Method first recovers the 
polynomial Q K (z) from the reflection response , r(n) , and then using reverse 
recursion , computes the reflection coefficients. But Qk(z) can be 
computed by the standard system identification problem also. 

The z-transform of r(n) , is given by , 



— P^zT" ’ ^cth the polynomials Q K (z) and P K (z) are polynomials of 
order equal to K. 

First , the AR-polynomial P«(z) <K=453> is computed from r(n) by 

solving the modified Yule-Walker equations given by , 
f r(K) r(K-i) .... rii) 1 


r(K+i) r<«) 

r(2) 


P(l) 


r(K+u 

r( 2 K— 1 ) r( 2 K— 2 ) .... 

ro<) 


PI2) 

as - 

r(K+2) 



P(K) 


r(2Ki 


where , p(n) are the autoregressive coefficients of the polynomial , P K (z). 
p(0)=i. The above modified equations are called extended , modified , or 

higher order Yule - Walker equations. These equations are solved by 
modified Levinson's recursive algorithm , which is used in the case of 
the matrix having Toeplitz structure. The number of multiplications needed 
for this modified Levinson's recursion is of the order of K 2 , but is 
about two times that required for the inversion of the Toeplitz matrix 
of order K. 

The MA - coefficients q(n) of the polynomial Q«(z) are computed by 
convolving p(n) with the reflection response , r(n>. The MA-coefficients 

oomputed from r(n) (obtained in step-3 of section 4.5.1) using this method 
is shown in Fig 4.5.4. The disadvantage of this method is that the 
order of the MA-polynomial , Q K (z) , must be known beforehand. If a 

wrong order is used , the computed MA-coefficients will not be correct. 

Also , this method gives only the MA-coefficients. Again the inverse 

recursion equations 3.43 are to be used to find the exact reflection 


coefficients . 



Fig. 4.4.1 

The real eepstrum of t2(r>), 
shown for n 200. 
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CHAPTER 5 

THE INFLUENCE OF RANDOM NOISE ON 
HOMOMORPHIC FILTERING 


Additive noise plays a critical role in homomorphic deconvolution to 
such an extent that the method becomes unreliable whenever the spectral 
amlitudes of the signal are very small over certain frequency bands and 
even a small amount of noise is present [83. That part of the noise 
which does not overlap the signal in the frequency spectrum can be 
eliminated by bandpass filtering. However the noise often overlaps the 
signal . 


The synthetic trace with additive noise is given by , 

x(n)s= t(n) + T?(n) (5.i> 

where, t(n) = s(n)#r(n)#p(n) and T?(n) is the additive noise. 

X(u) = T(u) + J?(w) 
log [ X(u) ] = log [ T(w) + 7?<u) ] 

= log [ T(u) ] + log [ i+ ] 


a log [ T(w) ] + 


where it is assumed that , 


Mu) 

T(u) 


« i 


(5.2) 

(5.3) 


Phase of X(u) ,the imaginary part of (5.2) is equal to , 


+ 


JM 

T(u) 


sin [ <Mw) - 3 ®.4) 

where, are the phases of the signal i.e. synthetic trace and the 

noise respectively . 


Thus it is clear from eqn,(5.4> that the phase of the trace 
becomes unpredictable and hence the phase unwrapping required in the 



computation of complex cepstrum will be unreliable due to the additive 
noise So the influence of random noise is quite critical in the 
computation of the complex cepstrum. In addition to this , noise is 
amplified during the computation of the cepstrum. The following factors 
influence the noise amplification C31 : 

i The nonlinear logarithmic operation exaggerates small components of the 
spectrum in the complex spectrum. 

2. Due to the fact that the noise is additive , it influences the complex 
cepstrum in a rather complicated way. 

3. The addition of random noise complicates primarily the phase spectrum , 

which has a strong influence on the complex cepstrum. However it is 
found that C31 the noise infuene is more severe on the high quefrency 

part of the cepstrum as compared to that around zero quefrencies. 

5.1 ILLUSTRATION : 

The influence of noise on homomorphic filtering is demonstrated for 

SNR =3, 10 and 40. 

1. White guassian noise is added to the synthetic trace t^n). The noisy 
trace for SNR = 40, 10 and 3 are shown in Figs. 5.1.1, 5.1.2 , 5.1.3 
respectively. 

2. The real and the complex cepstra for the exponentially weighted 

(a=0 .989) noisy trace, 5NR=40 are shown in Fig. 5.1.4, 5.15. From fig. 5.1.4 
it is clear that the complex cepstrum is influenced by the noise to a 
great extent even for high SNR (SNR=40). But we see that the low 

quefrency part of the real cepstrum is not very much affected by noise. 

3. Figs. 5.1.6, 5.1.7 & 5.1.8 show the recovered s m;ri (n) from the three noisy 



Figs. 5.1.6, 5.1.7 & 5.1.8 show the recovered s min (n) from the three noisy 


thaces 


It 

is 

observed 

that 

s(n) 

can not 

be 

recovered 

from 

the complex 

cepstrum 

for 

all the 

three 

cases, wheras 

the 

recovery 

of 

s(n) 

from 

real 

cepstrum 

is 

possible 

for 

all 

the three 

cases. However 

as 

the 

SNR 


decreases, the recovery also becomes poor. 
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Fig. 5.1.6 

Recovered emln(n) from the real 
cepstrum of the noisy tracefSNR— 40). 
dashed line - original emin(n). 






CHAPTER - 6 

CONCLUSION 


The main purpose of the present work is to try to solve the 
problems associated with the predictive deconvolution method. Synthetic 
traces are computed and then used to study the problems associated with 
the predictive deconvolution and to show the improvement in the results 
by combining predictive and homomorphic filtering. 

In random model , the predictive deconvolution method assumes s(n) 
to be minimum phase and r<n) to be an uncorrelated sequence. If sin) is 
a mixed phase sequence , then the recovered r(n) will contain an all-pass 
component. If r<n) is not a perfect white sequenoe , then its minimum 
phase equivalent sequenoe also gets deconvolved. To tackle these 

problems , homomorphic filtering is made use of along with the predictive 

filtering. The minimum phase equivalent of r(n) and the all-pass component 
of the mixed phase sin) are computed by homnomorphic filtering. Then the 
all-pass component is deconvolved from the recovered r(n) , followed by 

the restoration of the minimum phase equivalent of r(n). 


The 

success 

of the 

homomorphic 

filtering 

depends 

upon 

the 

exponential 

weighting 

required to 

make r(n) 

minimum 

phase. A 

trial 

and 


error method is proposed in this thesis. If the exponential weighting 
does not shift any of the zeros of s mix <n) lying outside the unit circle 
to inside the unit circle , then the maximum phase component of the 
exponentially weighted trace is the maximum phase component of s mix (n) 
itself. The all-pass component is computed from the maximum phase 



component of the s m ;y;(n^ This has an advantage that there is no need 

to choose the width of the cepstral window. On the other hand , if the 

exponentially weighting shifts some of the zeros of s mlx (n) to inside the 
unit circle , then the all-pass component corresponding to those zeros 
can not be computed by homomorphic filtering. 

In dynamic model , the dynamic predictive deconvolution (DPD) method 
can be applied to unit impulse response of the layered medium. So first 
homomorphic filtring is used to deconvolve s(n) from the trace. Then the 
DPD is used to recover the reflection coefficients of the layers. 

The effect of additive noise on Homomorphic deconvolution is 
studied. The detailed effects are not well understood C81, The noise tends 
to influence the long time ( i.e. high quefrency ) portion of the 

cepstrumC231 In this thesis this aspect has been observed in case of 
the noisy trace. Buttkus [31 has suggested that by high pass filtering 
the continous phase spectrum, the influence of the random noise at the 

deconvolution process can be made , reasonably small. He also mentions that 

though the experimental results confirms this, it is not known whether it 
is generally true. Thus , the problem , homomorphic deconvolution in 

presence of random noise, remains to be resolved C81. Further, better 

procedure for the determination of the exponential weighting needs to be 


evolved. 



APPENDIX 


HOMOMORPHIC SYSTEMS FOR CONVOLUTION 

Homomorphic systems are a class of nonlinear systems which satisfy 
a generalized principle of superposition. Any nonlinear system that 
satisfies the convolutional superposition property illustrated in Fig. A-i is 
called a homomorphic system for convolution. The procedure for 
homomorphic filtering is shown in Fig. A-2 and includes the following 
steps : 

til Computation of the complex spectrum of x(n) : 

The signal x(n) is the convolution of two signals , s(n) and r(n>. So 
the DFT of x(n) is , 

X<k) = S<k> R(k) = IX(k> l.e j6(k) . 

C23. LOG operation : 

The computation of the natural logarithm of the complex spectrum 
gives the additive superposition of the individual parts : 

logX(k)=logS(k)+logR(k) = log IX(k) I -f- 3 © ( k > 

[33. Computation of the complex cepstrum : 

The computation of the IDFT of the logarithm of the complex 

spectrum of x(n) gives the complex cepstrum of x(n). 

IDFTC logX(k>3 = x(n) = s(n> + rCn). 

Now the new "time domain" is called the ouefrency domain. The 

linear filtering of the complex cepstrum to suppress the undesired parts 

is called liftering. The reverse run through the the first three 
operations as shown in Fig. A-2 provides the final result of the 

homomorphic filtering in time domain. 



PROPERTIES OF COMPLEX CEPSTRUM : 

tl3 The complex cepstrum of a real sequence is also a real sequence. 

C23 The complex cepstrum decays at least as fast as i 

n - 

[33 If x(n) is of finite duration x<n) is nevertheless have infinite 

duration. 

[43 The complex cepstrum of a minimum phase sequence is zero for n<0 

[53. The complex cepstrum of a maximum phase sequence is zero for n>0. 

[63 The complex cepstrum of a minimum phase sequence satisfies 

x rnin^ n ^ ~ ( a ) Xmin(n) X m ;n(n-~ k) 

k=i 

[73. The complex cepstrum of a maximum phase sequence satisfies 

X n,ax<n> = *ma*<n> ~ £ ( x max (n; x max (n-k> 

k=n+i 

[83 The complex cepstrum of a pulse , whose spectrum is smooth tends to 

be concentrated around low quefrencies. 

[93. Let r(n) be a periodic sequence with period T >1 , defined by 

r(n)= 2 a k-S(n~kT) 
k=o 

where ty's are the impulse amlitudes. 

Then the complex cepstrum r(n) of r(n) is also a periodic impulse train 
with the same period T. 

C103. The removal of the first n complex cepstrum contributions of a 
sequence p(n) defined by 

p(n)»^<— a k S(n— kT) , 
k=o 

eliminates the first n time domain amplitudes of r<n) and reduces its 

remaining aplitudes to at most of their original amplitudes. 

n+i 

[113. Consider a minimum phase impulse train r<n> for which the duration 
between the first two nonzero samples is N , so that 

nj-riisN , n* arbitrary for k>2 



Then the 


t_ L-ifti.Ifey 


is zero for 0<n<N In 


cec’etf uir, of such a sequence 

general • rm* is nonzero only at quefrencies , 0 , n 2 _ ni , n s -n x , 

nm~nj , as well as at all positive linear combinations of these times. 

LU] Ulven a minimum delay sequence x(n) we form the reverse sequence 
x<-n) by folding x(n; about the origin. $<n) is the real oepstrum of x(n>. 

Then the real cepstrum of x(~n) = x<-n) 

[133 Given a minimum delay sequence b(n) we define the inverse sequence , 
b“ i( n) by 

b(n)*b -1 <n> = S(n). 

Then the real cepstrum of b“*(n) is equal to -b(n). 

[143 Real cepstrum of b -i (— n) is equal to -b(-n>. 

[153 A minimum phase sequence can be recovered from its real cepsrum 
(Fig. A 3). 

[163 The real cepstrum of the autocorrelation of a sequence is twice 

the real cepstrum of the sequence. 

[173 With x(n) expressed as the convolution of its minimum phase and 
maximum phase components , denoted as x miri (n) and x m a*(n) , respectively , 

x<n) = x min (n) + x max (n) 

[183. If xin) has a rational z-transform , then n.x(n) has a rational 
z-transform whose poles correspond to the poles and zeros of the z- 


transform of x(n). 
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